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FINAL  REPORT  N00014-01-1-0325:  Advancement  and  Application 
of  Multiphase  CFD  Modeling  to  High  Speed  Supercavitating  Flows 

Jules  W.  Lindau  Robert  F.  Kunz 


APPLICATIONS  AND  PAYOFF 

This  research  focuses  on  the  application  of  multiphase  CFD  modeling  for  high-speed 
underwater  flows.  The  payoff  for  this  and  other  related  investigations  will  be  analysis  input  for 
High  Speed  Supercavitating  Vehicle  (HSSV)  designs  including  test-bed  and  variable  cavitator 
designs.  The  tools  developed  will  be  used  for  assessment  of  candidate  cavitator,  fin  and  full-vehicle 
configurations.  They  will  also  provide  assessments  of  experiment  vs.  prototype  correspondence 
issues  such  as  integrated  passive  and  active  controllability,  water  tunnel  blockage,  and  scaling. 
Eventually  the  tools  will  be  evolved  to  where  they  can  reliably  provide  predictions  of  the  in  situ 
performance  capabilities  of  a  given  final  design,  assessment  of  the  limitations  of  a  given  design  for 
countermeasure  development  and  ultimately  analysis  of  maneuvering  vehicles. 

TECHNICAL  OBJECTIVES 

The  overall  technical  objective  of  this  program  has  been  to  develop  robust,  validated  CFD 
tools  for  supercavitating  flow  modeling  and  apply  the  analysis  capability  to  HSSV-relevant 
configurations. 


ABSTRACT 

Recent  progress  in  the  development  and  application  of  homogenous  multiphase  CFD 
methods  for  large-scale  gas  cavities  in  liquid  flows  are  presented.  The  focus  of  the  presentation  is 
on  work  in  n-species  transport  approaches  applied  to  developed  and  super-cavitation.  Numerical 
formulation,  physical  modeling,  and  applications  are  included.  Numerical  issues  to  be  discussed 
include:  preconditioning  in  the  context  of  incompressible  through  supersonic  Mach  numbers, 
arbitrary  numbers  of  species,  high  density  ratios,  and  large  gas  volume/mass  fractions; 
computational  grid  requirements;  dual-time  formulation;  and  general  large-scale  high-performance 
computing.  Physical  modeling  issues  to  be  discussed  include:  homogeneous  mixture 
compressibility;  mass-transfer  modeling;  condensable  and  non-condensable  gas  species;  turbulence 
modeling;  and  6-degree-of-freedom  [6DOF]  flowfield/rigid-body  interaction.  The  applications  to 
be  presented  range  from  naturally  cavitating  modeled  flow  on  simple  configurations  (ogives, 
nozzles,  airfoils/wedges)  to  more  industrially  relevant,  complex-geometry  applications  including 
turbomachinery  (cavitation  breakdown),  and  super-cavitation  (underwater  rockets,  hypervelocity 
darts,  condensable  and  non-condensable  cavities,  gas-on/off  transients,  twin-vortex  regime).  Recent 
applications  including  Detached  Eddy  Simulation  [DES]  of  natural  cavities,  and  6DOF  analysis  of 
a  propelled  notional  super-cavitating  vehicle  are  presented. 


1  INTRODUCTION 


Multiphase  flows  have  received  growing  research  attention  among  CFD  practitioners  due  in 
large  measure  to  the  evolving  maturity  of  single-phase  algorithms  that  have  been  adapted  to  the 
increased  complexity  of  multicomponent  systems.  However,  there  remain  a  number  of  numerical 
and  physical  modeling  challenges  that  arise  in  multi-phase  CFD  analysis  beyond  those  present  in 
single-phase  methods.  Principal  among  these  are  large  constituent  density  ratios,  the  presence  of 
discrete  interfaces,  significant  mass  transfer  rates,  non-equilibrium  interfacial  dynamics, 
compressibility  effects  associated  with  the  very  low  mixture  sound  speeds  which  can  arise,  the 
presence  of  multiple  constituents  (viz.  more  than  two)  and  void  wave  propagation.  These  naturally 
deserve  special  attention  when  a  numerical  method  is  constructed  or  adapted  for  multi-phase  flows. 

The  ability  to  properly  model  multiphase  flows  has  significant  potential  engineering  benefit.  In 
particular,  liquid-gas  flows  where  a  persistent  gas  phase  occupies  a  significant,  large-scale,  volume 
are  considered.  Large-scale  vaporous  cavities  may  appear  in  flows  about  submerged  objects  due  to 
reduction  in  local  static  pressure  below  the  local  vapor  pressure.  Such  flows  represent  traditional 
cavitation.  In  addition,  ventilated,  gaseous  cavities  are  present  in  specialized  high  speed 
applications  where  enhanced  cavities  are  used  to  improve  hydrodynamic  performance.  A  third  type 
of  cavitating  flow  occurs  when  liquid  and  gas  temperatures  are  elevated  significantly  due  to 
stagnation  conditions,  chemistry,  or  other  heat  sources.  In  this  case  the  equilibrium  vapor  pressure 
may  be  affected  by  the  local  temperature.  Traditionally,  cavitation  has  had  negative  implications 
associated  with  damage,  reduction  in  hydrodynamic  performance,  and/  or  noise.  However,  for 
certain  applications,  a  natural  or  ventilated  supercavity  has  demonstrated  potential  benefit. 
Intentionally-placed  supercavities  have  long  been  applied  in  the  design  of  hydrofoils  as  well  as  high 
performance  marine  propellers.  An  application  of  particular  interest  to  the  authors  is  that  of  a  very 
high  speed  aquatic  vehicle  mostly  enshrouded  in  a  supercavity,  and,  thus,  experiencing  a 
significantly  reduced  drag  relative  to  a  fully-wetted,  counterpart. 

Underwater  multiphase  flows  typically  involve  the  presence  of  compressibility  effects  in  a  flow 
that  is  largely  incompressible.  For  typical  marine  applications  this  is  due  to  the  well-known 
dramatic  decrease  of  the  acoustic  speed,  from  either  isolated  phase,  for  a  homogeneous,  liquid-gas 
mixture  [1].  The  situation  is  illustrated  in  Fig.  1  for  a  saturated  water  mixture  at  300K.  The 
homogeneous  mixture  relation  for  isothermal  sound  speed  has  been  plotted  as  a  function  of  vapor 
volume  fraction.  Note  that,  for  moderate  vapor  volume  fractions,  the  sonic  velocity  is  less  than 
5m/s,  which  is  almost  two  orders  of  magnitude  less  than  the  sonic  velocities  in  either  of  the  pure 
phases.  As  has  been  noted  in  many  other  sources,  in  the  presence  of  cavitation,  velocities  normally 
associated  with  subsonic,  incompressible  flows  may  result  in  supersonic  conditions  and  associated 
strong  wave  formations  such  as  shocks.  It  is  therefore  evident  that  such  compressible  phenomena 
must  be  represented  accurately  in  multiphase  simulations. 

The  class  of  multiphase  flows  under  consideration  here  is  developed  and  supercavitating  flows, 
wherein  significant  regions  of  the  flow  are  occupied  by  gas  phase,  thus  the  terminology,  “large- 
scale.”  Depending  on  the  configuration,  such  “developed”  [2]  cavities  are  composed  of  a  mixture 
vapor  and  non-condensable  gases.  Historically,  most  efforts  to  model  large  cavities  relied  on 
potential  flow  methods  applied  to  the  liquid  flow,  while  the  bubble  shape  and  closure  conditions 
were  specified.  Adaptations  of  potential  flow  methods  remain  in  widespread  use  today  [3],  due  to 
their  inherent  computational  efficiency,  and  their  proven  effectiveness  in  predicting  numerous  first 
order  dynamics  of  super-cavitating  configurations.  Recently,  more  general  CFD  approaches  have 
been  developed  to  analyze  these  flows.  In  one  class  of  methods,  a  single  continuity  equation  is 
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considered  with  the  density  varying  abruptly  between  vapor  and  liquid  densities  through  an 
equation  of  state.  Such  “single-continuity-equation-homogeneous”  methods  have  become  fairly 
widely  used  for  sheet  and  super-cavitation  analysis  ([4-8],  for  example).  Although  these  methods 
can  directly  model  viscous  effects,  they  are  inherently  unable  to  distinguish  between  condensable 
vapor  and  noncondensable  gas,  a  requirement  of  ventilated  super-cavitating  vehicle  analysis.  By 
solving  separate  continuity  equations  for  liquid  and  gas  phase  fields,  one  can  account  for  and  model 
the  separate  dynamics  and  thermodynamics  of  the  liquid,  condensable  vapor,  and  noncondensable 
gas  fields.  Such  multispecies  methods  are  also  termed  homogeneous  because  interfacial  dynamics 
are  neglected,  that  is,  there  is  assumed  to  be  no-slip  between  constituents  residing  in  the  same 
control  volume.  A  number  of  researchers  have  adopted  this  level  of  differential  modeling,  mostly 
for  the  analysis  of  natural  cavitation  where  two  phases/constituents  are  accounted  for  ([9-1 1],  for 
example).  This  is  the  level  of  modeling  employed  here,  though  a  multiple  species  formulation  is 
used  to  account  for  N-gaseous  fields.  For  two  phases/constituents  these  methods  are  very  closely 
related  to  the  “single-continuity-equation-homogeneous”  methods  addressed  above  with  interfacial 
mass  transfer  modeling  supplanting  an  equation  of  state. 

Full-two-fluid  modeling,  wherein  separate  momentum  (and  in  principle  energy)  equations  are 
employed  for  the  liquid  and  vapor  constituents,  have  also  been  used  to  model  natural  cavitation 
[12].  However,  in  sheet-cavity  flows,  the  gas-liquid  interface  is  known  to  be  nearly  in  dynamic 
equilibrium;  for  this  reason,  we  do  not  pursue  a  full  two-fluid  level  of  modeling. 

Developed  cavitating  flows  are  characterized  by  large  density  ratios  (liquid  to  vapor  density 
ratios  of  greater  than  104  is  typical  of  marine  cavitation),  relatively  discrete  cavity-free  stream 
interfaces  and,  due  to  ventilation,  multiple  gas  phase  constituents.  Accordingly,  the  CFD  method 
employed  must  accommodate  these  physics  effectively.  Most  relevant  applications  exhibit  large 
scale  unsteadiness  associated  with  re-entrant  jets,  periodic  ejection  of  non-condensable  gas,  cavity 
“pulsations”,  and  large  scale  vortices.  In  addition,  we  are  interested  in  capturing  the  flow  about 
maneuvering  vehicles  with  articulated  control  surfaces.  Accordingly,  we  and  others  ([7, 10, 1 1, 13]) 
employ  a  time-accurate  formulation  in  the  analysis  of  large  scale  cavitation. 

In  the  presentation  that  follows,  the  governing  multiphase  physical  system  will  be  defined.  A 
local  preconditioning  methodology,  essential  to  rendering  model  solutions  efficiently  and 
accurately,  will  be  reviewed.  The  modification  to  turbulence  transport  modeling  known  as 
Detached  Eddy  Simulation  (DES),  as  applied  to  multiphase  model  solutions  presented  here,  will  be 
described.  The  method  of  coupling  rigid  body  motion  to  multiphase  fluid  flow  as  applied  to  model 
solutions  will  also  be  described.  The  methods  of  numerical  integration  of  the  governing  equations 
will  be  given,  and  finally,  a  section  of  computational  results  is  presented.  Results  supporting  the 
validity  and  documented  the  capabilities  of  the  modeling  method  have  been  chosen. 


NOMENCLATURE 

Ci,  C2 

turbulence  model  constants 

cq 

dimensionless  ventilation  rate  [q/(UoDc2)l  (q  is  volume  flow  rate) 

Cprod,  Cdest 

mass  transfer  model  parameter 

cD 

specific  heat  and  pressure  coefficient 

c,  c' 

sound  speed,  pseudo-sound  speed 

Dc 

cavitator  diameter 

Froude  Number  [Uo/(Lg)1/2l 

e 

h,  h0 
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k,  s 

turbulence  model  eddy  viscosity  and  dissipation  rate 

L. 

turbulent  length  scale 

m+,m~,d> 

nass  transfer  source  terms 

P 

turbulence  model  production  term 

Pr 

Prandtl  Number 

&k->  CTe 

turbulence  model  Schmidt  Numbers  for  k  and  e 

P,q,r 

vehicle  relative  rotation  rates  (also  p  is  thermodynamic  pressure) 

Pv 

saturation  vapor  pressure 

T 

temperature 

S 

rigid-body  motion  dependent  variable  vector 

too 

time  parameter  in  mass  transfer  model 

M/',  LIco 

Cartesian  velocity,  freestream  velocity 

Xi 

Cartesian  coordinate 

Y,  7 

mass  fraction,  relative  (to  total  gas)  mass  fraction 

+ 

y 

dimensionless  wall  distance  (boundary  layer) 

CX/,Ot  £?(Xy 

liquid,  gas,  and  vapor  volume  fractions 

rp 

preconditioning  matrix 

Ye 

specific  gravity  of  vehicle 

8// 

Kronecker  delta 

f^m->  ft  m.t 

mixture  molecular,  and  eddy  viscosity 

K  m.t 

mixture,  turbulent  conductivity 

p,p„pg,pv 

mixture,  and  isolated  liquid,  gas,  and  vapor  densities 

cr 

cavitation  index  [(poo-pv)/(0.5p/UoO2)] 

2  MODEL  FORMULATION 

The  physical  model  is  locally  homogenous,  multifield  and  single-fluid.  Thus  the  model  consists 
of  species  conservation  equations  coupled  to  a  single  momentum  and  energy  equations.  The  form  is 
further  refined  by  a  time-dependent  (unsteady),  Reynolds-Averaged  Navier-Stokes  (URANS) 
model  using  a  two-equation  turbulence  closure.  For  some  applications  turbulence  modeling  is 
modified  via  a  DES  technique  [14].  Mass  transfer  modeling  between  liquid  and  vapor  phases  is 
achieved  with  an  ad-hoc  model  described  in  this  section.  Although  no  results  including  gas-gas 
mass  transfer  are  presented.  Gas  phase  chemistry  modeling,  based  on  the  work  of  Venkateswaran 
and  others  [15],  follows  methodology  similar  to  the  liquid-vapor  model  presented  here.  Methods  of 
averaging  the  mixture  scalar  quantities  are  based  on  mass  for  intensive  quantities  and  volume  for 
all  others.  Quantities  of  interest  are  specific  heat,  molecular  and  eddy  viscosity,  laminar  and 
turbulent  Prandtl  number,  and  conductivity. 

The  differential  form  of  the  computational  model  in  Cartesian  tensor  notation  is  given  in  Eqn.  1. 
The  corresponding  conservative  variables,  primitive  variables,  flux  vectors,  and  source  terms  are 
defined  in  Eqn.  2.  This  equation  differs  in  appearance  from  the  basic  physical  system  due  to 
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appearance  of  the  term,  Tp  — ,  a  result  of  the  Choi-Merkle,  dual-time  preconditioned  approach 

dr 

[15].  Definition  of  the  preconditioning  matrix  is  given  in  Subsection  2.1,  a  description  of  the  DES 
method,  as  applied  here,  is  given  in  Subsection  2.2,  and  a  description  of  the  6-degree-of-freedom 
coupled  fluid-rigid-body-motion  methodology  is  given  in  Subsection  2.3. 

d&+Y„dQ  +  F  _Fv  =H  (1) 

8t  8t  JJ  jJ 


Q= 


~  p~ 

a,p, 

a,p,Uj 

UJ 

pUj 

P",«j  +  Psij 

T 

Qc  = 

e 

FJ= 

Uj(e  +  p) 

FJ  = 

ag 

pg^X 

Psaj'uj 

kJ 

PgaX. 

pgagY*Uj 

0 


»itij  +  KmJ 

0 
0 


dT 

dx. 


m+  +  m 

o, 

H  = 

0 

-(m+  +m~) 

(2) 


In  the  vector  form  of  Eqn.  2,  the  first  equation  represents  continuity  of  the  liquid-phase  written 
in  terms  of  volume  fraction.  The  next  two  equations  are  momentum  and  energy  equations  for  the 
mixture.  The  last  two  equations  deserve  further  comment.  They  represent  individual  gaseous 
species  continuity  equations.  The  subscript  g  is  used  to  denote  all  gaseous  (including  vapor) 
species,  while  the  subscript  v  denotes  the  isolated  vapor  phase  and  k  denotes  non-condensable 
gaseous  species.  For  generality,  k=l,...,N,  where  N  represents  the  number  of  non-condensable 
gaseous  species.  The  tilde  above  the  density  indicates  that  the  quantity  is  defined  by  the  rule  of 
Amagat  rather  than  Dalton,  i.e.  considering  the  isolated  constituent,  rather  than  the  mixture.  Note 
that,  the  mass  fraction  and  the  volume  fraction  are  related:  pYg  -  ccgpg .  Further,  the  mixture 

density  is  given  by  Eqn.  3. 

N  (3) 

P  =  agpg  +  a>Pi  =  0  -  ai)pg  +  a,p,  =  +  IX  Pk  +  aiPi 


k= 1 


The  superscript,  *,  indicates  a  relative  gas  mass  fraction,  defined  in  Eqn.  4: 


K*  =  —  =  - 


and  T  =  —  = 

T„ 


(4) 


j;+IX  /g  n+l X 

k= 1  *=1 

The  viscous  stress  tensor  takes  on  the  usual  form. 

hj  =  Mn,.,  [u>,j  +  » j.i  ~  Sij  j 


(5) 


In  Eqn.  6,  necessary  mixture  quantities  are  defined.  The  subscript  m  indicates  a  liquid-gas 
mixture  quantity.  The  mixture  viscosity  and  Prandtl  number  are  computed  based  on  a  local  volume 
average.  The  mixture  specific  heat,  is  computed  based  on  a  mass  average. 
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Finally,  two-equation  turbulence  closure  is  employed. 
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The  transfer  of  mass  between  liquid  and  vapor  states  is  handled  with  simple  finite-rate  relations 
given  in  Eqn.  8.  The  mass  transfer  rate  between  gaseous  species  due  to  chemical  reactions  is 
indicated  in  Eqn.  1  by  a>k . 
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The  mass-transfer  model  is  presented  without  rigor.  It  is  supported  by  its  supposed  relation  to 
the  understood  physics:  Destruction  of  liquid  is  related  to  the  difference  of  the  local  pressure,  p,  and 
vapor  saturation  pressure,  pv,  and  the  production  of  liquid  is  based  on  an  interpretation  of  the  work 
of  Hohenberg  [16].  Further  support  is  afforded  by  the  presented  computational  results.  Within  the 
context  of  Eqn.  8,  saturation  vapor  pressure  is  an  assumed  function  of  temperature;  thermal  effects 
on  phase  change  may  be  included.  Cdest  and  Cpro d  are  empirical  constants.  Unless  otherwise  noted, 
for  all  results  presented,  Cdest=Cprod=100,  and  the  time  scale  parameter,  tx,  is  based  on  a  convective 
time  scale  roughly  related  to  cavity  formation.  For  modeled  cases,  where  a  cavitator  is  well  defined, 
this  time  scale  has  been  based  on  the  cavitator  diameter  and  the  free  stream  velocity.  For 
turbomachinery  and  hydrofoil  model  problems,  this  was  based  on  the  maximum  chord  length  and 
the  relative  inflow  velocity.  As  has  been  presented  in  Ref.  [17],  model  results  are  relatively 
insensitive  to  these  choices. 

In  the  presented  work  all  walls  are  adiabatic.  The  model  system  is  closed  with  the  inclusion  of 
equations  of  state  for  the  constituent  phases  and  species.  For  the  liquid  phase,  for  instance, 
p,  -  p,(p,T)  and  h,  =h,(p,T).  For  the  constituent  gaseous  species,  including  the  vapor,  the 

species  density  and  enthalpy  may  be  expressed  similarly.  Typically  an  ideal  gas  equation  of  state  is 
employed. 
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2.1  Local  Preconditioning 

In  order  to  obtain  convergent  and  accurate  numerical  solution  of  a  linear  algebraic  system,  the 
condition  number  of  that  system  must  be  sufficiently  small.  An  appropriate  definition  of  the 
condition  number  is  the  magnitude  of  the  ratio  of  the  largest  magnitude  eigenvalue  to  the  smallest 
magnitude  eigenvalue  of  the  matrix  system  [18].  The  optimal  value,  therefore,  is  one.  Using  the 
current  finite-volume  method  of  discretization,  the  global  algebraic  system  may  be  divided  into 
local  subsystems  considering  the  solution  of  equations  related  to  individual  control  volumes.  In 
order  to  obtain  a  well-conditioned  system  of  equations  for  solution  of  multiphase,  turbulent, 
viscous,  compressible  flow  problems  of  interest,  it  is  beneficial  to  examine  and,  by  rescaling  the 
time-derivative  terms,  manipulate  the  inviscid  eigensystem.  The  condition  of  this  subsystem  tends 
to  indicate  the  condition  status  of  the  more  complex  full  numerical  system.  This  has  been  explained 
in  detail  in  Ref.  [15]  and  is  the  procedure  followed  here.  The  one-dimensional  form  of  the 
governing  multi-phase  equations  is  considered.  Extension  to  multiple  dimensions  is  straight¬ 
forward  and  does  not  effect  the  conclusions  drawn  here. 

da+3£=0  (9) 

dr  dx 

A  set  of  primitive  variables  may  be  chosen  with  the  appropriate  transformation  Jacobian.  Here 
the  transformation  is  not  exact  but  is  a  preconditioning  matrix  containing  pseudo-properties 
substituted  for  the  physically  correct  definitions.  These  substitutions  will  be  shown  to  favorably 
affect  accuracy  and  convergence  of  the  solution  method  without  altering  the  correct  physics  of  the 
solution.  This  form  of  the  equations  is  given  in  Eqn.  10.  In  Eqn.  1 1,  the  conservative  variables,  the 
fluxes,  and  two  primitive  variable  sets  are  defined.  One  is  denoted  by  the  subscript  2  and  is 
convenient  for  development  of  the  preconditioning  matrix.  The  other  is  indicated  by  subscript  1  and 
is  preferred  (by  the  authors)  for  computations.  The  elements  of  these  sets  are  chosen  for  algebraic 
convenience  in  developing  the  preconditioning  matrix. 

rpd&+dE=Q  (10) 

2  dr  dx 

pyi~\  r  pyiu 

pii  pu 2  +  p 

Qc=  e  E=  u(e  +  p )  Q2 

PYv  pYvu 

_pYk\  [  PYk“  . 

The  preconditioning  matrix  is  defined  in  Eqn.  12.  This  formulation  has  been  given  in  Ref.  [19]. 
A  similar  but  more  detailed  preconditioning  derivation  for  the  2-phase,  liquid-vapor  system,  was 
given  in  Ref.  [20]. 
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The  corresponding  flux  Jacobian,  A2  = - ,  is  given  in  Eqn.  13. 
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In  each  of  these  matrices,  in  the  last  row  and  last  column,  indexing  over  the  number  of 
noncondensable  gaseous  species  is  implied.  Thus  the  definition  extends  to  an  arbitrary  number  of 
species.  The  eigenvalues  of  the  preconditioned  system  are  given  in  Eqn.  14. 
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The  preconditioned,  pseudo-sound  speed,  c' ,  is  given  in  Eqn.  15.  In  Eqn.  14,  c  is  the  physical 
sound  speed.  Note  that,  c,  the  isentropic,  frozen  mixture  sound  speed  is  given  by  removing  the 
primes  from  Eqn.  15. 
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An  appropriate  preconditioner  definition  for  the  pseudo-sound  speed  is  the  magnitude  of  the 
local  convective  velocity.  The  appropriateness  of  this  definition  may  be  confirmed  by  applying 
Eqn.  14.  The  resulting  eigensystem  is  well-conditioned  for  all  Mach  numbers  independent  of 
important  effects  such  as  phasic  density  ratio  and  sound  speed.  As  is  the  case  with  low-Mach 
number  preconditioning  for  single-phase  flows,  in  numerical  practice,  the  pseudo-sound  speed 
value  should  be  lower  limited  [15].  For  applications  presented  here,  this  lower  limit  was  chosen 
based  on  the  free  stream  velocity  magnitude.  For  cases  where  the  free-stream  phase  is  modeled  as 
incompressible,  the  optimal  lower  limit  was  found  to  be  the  square-root  of  ten  times  the  free  stream 
magnitude.  For  all  other  cases,  the  lower  limit  was  set  equal  to  that  magnitude.  Other  more  rigorous 
methods  of  limiting  the  local  preconditioning  definition  have  been  derived  in  Ref.  [21].  These 
definitions  have  been  applied  with  success  to  the  physical  system  outlined  here.  However,  no 
results  based  on  this  methodology  are  presented.  Another  requirement  is  that  the  preconditioned 
system  should  revert  to  the  physically  exact  system  when  the  local  velocity  magnitude  is 
supersonic.  In  this  case  the  physical  system  is  automatically  well  conditioned.  Therefore,  in 
computational  practice,  the  physical  sound  speed  is  substituted  for  the  pseudo-sound  when  the  local 
Mach  number  is  greater  than  one.  This  should  yield  a  continuous  transition  from  preconditioned  to 
physical  system  at  the  sonic  point.  Further  discussion  of  the  appropriate  choice  of  pseudo¬ 
properties  is  given  in  Ref.  [15], 

The  computational  form  chosen  is  given  by  primitive  variable  set  1 .  Therefore,  transformation  is 
necessary.  The  transformation  is  obtained  using  the  Jacobian  relating  variable  set  1  to  set  2  as 
shown  in  Eqns.  16  and  17. 

iySfc+*-o  <l6) 

dr  dx 


T^p  _  ^Ql  T'P 

1  ~  a  n  2 

°Q\ 

The  transformation  Jacobian  is  given  by  Eqn.  1 8. 
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2.2  Detached  Eddy  Simulation 


Detached  Eddy  Simulation  (DES)  and  other  hybrid  RANS-Large  Eddy  Simulation  (LES) 
approaches  have  emerged  recently  as  a  potential  compromise  between  RANS  based  turbulence 
models  and  full  LES.  Specifically,  for  flows  characterized  by  large  scale  separation  and 
unsteadiness,  DES  has  been  shown  by  a  number  of  workers  to  date  to  provide  improved  mean  flow 
predictions  [14, 22, 23]. 

DES  is  a  generalization  of  RANS  modeling.  The  underlying  turbulence  model  for  DES  is  a 
RANS  subgrid  model  with  the  characteristic  length  scale  reinterpreted  in  terms  of  the  local  grid 
scale.  When  the  local  mesh  is  fine  relative  to  the  turbulence  mixing  length,  the  DES  model 
becomes  an  LES  with  a  Smagorinsky-like  subgrid  closure  [14].  In  these  resolved  regions,  modeling 
error  is  minimized.  The  solution  remains  RANS  elsewhere  and  retains  the  solution  efficiencies  of 
RANS. 


DES  is  formulated  for  and  is  particularly  suited  to  regions  of  massive  separation  where  RANS 
closures  are  challenged.  With  minimal  increases  in  grid  density,  the  governing  equations  are 
allowed  to  support  the  most  energetic  turbulence  scales  yielding  improved  solution  fidelity.  Near 
boundaries,  where  kinematic  constraints  limit  the  scale  of  turbulence  eddies  and  where  grid 
resolution  demands  for  traditional  LES  are  most  demanding,  DES  reverts  back  to  RANS.  By 
design,  RANS  models  can  perform  very  well  near  solid  surfaces. 

In  the  context  of  the  k-e  model  given  in  Eqn.  7,  a  DES  variant  is  obtained  by  defining  a  modified 
turbulence  energy  destruction  term  as  (see  [14]  for  example): 

-Pe  ~Pe  ■  fdes  (19) 
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2.3  Fluid-Rigid-Body  Coupled  Motion 

The  multiphase  fluid  dynamics  model  is  coupled  to  the  general  equations  of  vehicle/rigid  body 
motion  (6DOF).  Thus  the  unsteady  flow  (URANS)  and  the  6DOF  are  submodels  to  be  coupled  for 
solution  of  the  overall  model  problem.  Coupling  between  the  submodels  is  achieved  via  modeling 
the  shear  and  normal  forces  on  the  surface  of  the  rigid  body.  The  shear  and  normal  forces  are 
results  of  the  instantaneous  solution  to  the  URANS.  In  this  model,  the  normal  forces  include  both 
the  pressure  on  solid  surfaces  and  the  impulse  through  propulsion  and  ventilation  computational 
boundaries.  Shear  forces  are  computed  based  on  a  turbulent  wall  function.  The  URANS  forces  are 
then  resolved  into  appropriate  instantaneous  force  and  moment  components  for  the  6DOF.  The 
result  of  the  6DOF  solution  is  then  translation  and  rotation  of  the  rigid  body.  With  adequate 
discretization,  the  6DOF  and  the  URANS  then  tend  to  the  proper  physics.  The  current  presentation 
and  results  are  for  a  single,  completely  rigid,  body  and  closely  follow  the  methodology  of  Dreyer 
[24]  for  single-phase  flow.  Formulation  of  the  single-phase  pseudocompressible  formulation  for 
nonstationary  grid  systems  has  been  well  presented  by  Taylor  [25].  Solutions  with  multiple  bodies, 
or  a  single  body  with  articulated  control  surfaces,  represent  straightforward  extensions  of  the 
current  methodology. 

Computation  of  fully  coupled  multiphase  URANS-6DOF  modeling  is  achieved  with  a 
computational  flow  grid  fixed  to  the  rigid  body.  Flow  solution  is  then  obtained  in  an  inertial 
reference  frame.  The  rigid  body  motion  problem  is  solved  in  a  body-fixed  reference  frame.  These 
choices  follow  the  work  of  Dreyer  et  al.  [24].  During  each  physical  time  step,  the  hydrodynamic 
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forces  on  the  body  are  computed  based  on  the  instantaneous  flow  field  solution,  then  the  position  of 
the  body,  is  updated  based  upon  solution  of  the  6DOF  equations,  and  then  the  flow  solution  is 
updated.  More  detail  related  to  this  method  is  given  in  Ref.  [26]. 

2.4  Numerical  Integration 

The  described  model  flow  and  rigid  body  motion  equations  are  solved  in  the  UNCLE-M  code. 
The  flow  solver  is  structured,  multiblock,  implicit  and  parallel  with  upwind  flux-difference  splitting 
for  the  spatial  discretization  and  Gauss-Seidel  relaxation  for  the  inversion  of  the  implicit  operator. 
Primitive  variable  (MUSCL)  interpolation  with  limiting  based  on  the  solution  volume  fraction  is 
applied  to  retain  higher  order  accuracy  in  flow  fields  containing  physical  discontinuities.  Source 
term  linearization  is  approached  on  a  case-by-case  basis.  However,  it  has  been  found  that  for  flows 
with  free-stream  phase  Mach  numbers  approaching  the  incompressible  limit,  in  keeping  with 
previous  work  [27],  those  source  terms  associated  with  vapor  production  should  be  linearized  for 
inclusion  in  the  implicit  linear  system  left-hand-side.  Terms  associated  with  liquid  production  are 
treated  explicitly  and  under-relaxed  with  a  factor  of  0.1.  For  cases  with  compressible  free-stream 
Mach  numbers,  it  has  been  found  that,  depending  on  initial  conditions,  full  source-term 
linearization  is  appropriate.  With  typical  initial  conditions,  until  the  residual  drops  one  to  two 
orders  of  magnitude,  under-relaxation  with  no  source-term  linearization,  combined  with  reduced 
integration  step  size,  may  be  the  best  approach.  After  this  start-up  effort,  full-linearization  with 
larger  integration  steps  usually  is  effective.  This  procedure  is  in  accordance  with  the  work  of 
Venkateswaran  et  al.  [15].  At  each  pseudo-time  step,  the  turbulence  transport  equations  are  solved 
subsequent  to  solution  of  the  mean  flow  equations.  Further  details  regarding  the  numerical  method 
are  available  in  Ref.  [27].  For  cases  with  rigid  body  motion,  the  6DOF  model  is  a  set  of  six 
ordinary  differential  equations  in  time  and  is  integrated  with  a  five-stage,  fourth-order  Runge-Kutta 
method.  The  same  integration  method  was  used  by  Dreyer  et  al.  [24]. 

3  APPLICATIONS 

Several  multiphase,  compressible  and  incompressible,  two  and  three-phase  flowfields  have  been 
modeled.  These  flows  include  validative  cases  compared  to  experimental  results,  as  well  as 
demonstrative  cases  highlighting  physically  sensible  results  and  capabilities  of  the  modeling 
method. 

3.1  Cavitating  Flow  in  a  Venturi 

Reboud,  Stutz,  et  al.  [28]  have  performed  detailed  unsteady,  flowfield  measurements  of 
vaporous  cavitating  flow  in  the  two-dimensional  Venturi  section  of  a  water  tunnel.  The  test  section 
captures  significant  physics  found  on  the  suction  side  in  a  blade  passage.  Thus  their  experiment  and 
the  current  model  results  represent  partially  cavitating  flow  in  a  turbomachinery-like  environment. 
In  Fig.  2,  the  average  and  RMS  fluctuating  portions  of  the  liquid  volume  fraction  is  presented  based 
on  the  modeled  flow.  This  figure  serves  to  illustrate  the  geometry  of  the  modeled  test  section  as 
well  as  the  results  obtained  during  modeling.  The  test  section  had  a  height  at  the  throat  equal  to 
43.7  mm  and  a  constant  width  equal  to  44  mm.  The  nominal  cavity  length  for  comparison  here  was 
80  mm  in  the  horizontal  direction.  The  nominal  angle  of  the  lower  surface  of  the  Venturi 
downstream  of  the  throat,  with  respect  to  a  horizontal  line,  was  4°.  The  experiments  were 
conducted  at  Reynolds  numbers  based  on  cavity  length  from  4.3  xlO5  to  2.1><106  and  at  a  range  of 
cavitation  numbers,  based  on  the  upstream  pressure  and  velocity,  from  0.6  to  0.75.  It  may  be  seen 
from  the  figure  that  although  there  is  a  high  degree  of  unsteadiness  in  the  region  of  the  cavitating 
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flow,  this  unsteadiness  is  confined  to  the  test  section  area.  This  is  consistent  with  the  commentary 
[28]. 

In  Fig.  3,  the  current  computational  results  are  presented  with  the  experimental  data  [28].  In  the 
figure,  results  have  been  plotted  at  the  five  measurement  stations  used  in  the  experiments.  Each  of 
these  stations  is  given  at  a  horizontal  position.  The  experimental  cavity  was  initiated  due  to  the 
suction  peak  on  the  lower  surface  of  the  throat  of  the  test  section,  a  reference  position  of  x=0.  In 
part  (a)  of  the  figure,  the  mean  vapor  volume  fraction  is  plotted  at  the  five  axial  stations.  Clearly  the 
model  tends  to  over  estimate  the  void  fraction,  particularly  at  the  forward  region  of  the  cavity,  at 
x=22.5mm.  However,  the  average  quantities  are  in  excellent  agreement  at  the  tail  end.  Similar 
agreement  is  demonstrated  with  the  unsteady  portion  of  the  RMS  void  fraction,  part  (b).  Here,  the 
error  is  greater  in  the  closure  region,  at  x=60mm  and  x=80mm.  Considering  the  difficulty  of 
modeling  in  the  closure  region,  this  level  of  agreement  is  also  pleasing.  In  part  (c),  the  average  axial 
velocity  is  given,  and,  in  part  (d),  the  RMS  fluctuating  component  is  given  at  the  five  measurement 
stations.  Here  the  agreement  is  generally  good,  except  at  the  tail  end  for  the  average  velocity,  where 
the  reverse  flow  is  missed.  It  should  be  noted  that,  by  application  of  a  two-phase  Navier-Stokes 
model  based  on  a  barotropic  state  law,  Reboud  et  al.  [28]  were  able  to  obtain  similarly  good 
agreement  with  the  experimental  data. 

3.2  Naturally  Cavitating  Flow  over  Axisymmetric  Bodies 

Turbulent,  naturally  cavitating  flow  over  axisymmetric  bodies  is  known  to  be  a  highly  nonlinear 
and  three-dimensional  event.  This  is  clearly  illustrated  in  Fig.  4.  Here,  a  photograph  during  water 
tunnel  testing  of  a  blunt  cavitator  at  zero  angle-of-attack,  o~0.35,  and  /?ez>»1.5xl05  is  shown  in  part 
(a)  to  be  compared  and  contrasted  to  the  model  result  in  part  (b).  To  obtain  the  model  result, 
turbulent  vaporous  cavitating  flow  over  a  blunt  cavitator  was  modeled,  a  was  set  to  0.4  and  Reo 
was  1.46xl05.  An  appropriate  high  Reynolds  number  grid  with  approximately  1.2  million  nodes 
was  used.  The  snapshot  of  part  (b)  represents  a  physical  time  slice  taken  after  a  clear  model  cavity 
cycle  had  been  established.  In  the  figure,  an  isosurface  at  ap0.5  has  been  presented  with  selected 
streamlines,  and  the  surface  of  the  cylinder  has  been  colored  by  volume  fraction.  The  streamlines 
are  merely  suggestive  (but  helpful),  as  they  have  been  generated  based  on  instantaneous  velocity 
vectors.  Clearly  in  neither  the  model  result  nor  the  photograph  is  the  flowfield  in  and  around  the 
cavity  axisymmetric.  It  is  suspected  that  physical,  chaotic,  dynamic  interdependencies  are 
responsible.  Via  an  adequate  level  of  modeling,  a  characteristic  of  the  real  flow  has  been  well 
captured.  Positive  understanding  of  the  causal  mechanisms  is  a  subject  for  further  research. 

Dependence  of  the  nature  of  the  result  on  the  method  of  numerical  integration  has  been 
determined  to  be  minimal.  This  was  done  by  numerical  experiment.  The  object  of  this  experiment 
was  to  determine  if  the  helical  reentrant  jet  procession  was  dependent  on  the  current  method  of 
integration,  even  though  the  computational  grid  was  axially  symmetric.  A  portion  of  the  result  is 
shown  in  Fig.  4c.  Possible  bias  was  hypothesized  to  come  from  such  factors  as  initial  Gauss-Seidel 
sweep  direction,  etc.  The  entire  model  state  from  the  unsteady  integration  results  for  a  steadily 
clockwise  processing  jet  was  transformed  to  a  mirror  image  of  that  state.  This  physical  mirror  state 
was  then  used  as  initial  conditions  to  complete  the  experiment.  All  other  parameters  were  unaltered, 
including  the  structured  grid  numbering  and  computational  (ij,k)  coordinates.  Integration  results  on 
both  the  mirror  and  the  original  solutions  were  then  compared.  The  solutions  were  found  to  be 
mirror  images  of  each  other.  Solutions  had  statistically  identical  but  opposite  helical,  reentrant 
flows;  one  processed  clockwise  and  the  other  in  the  counterclockwise  direction.  All  other 
quantitative  comparisons  showed  identical  but  mirror-image  solutions. 
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3.3  Propeller  Cavitation  Breakdown  Analysis 

A  no-skew  unducted  propeller,  designated  P4381,  was  tested  by  Boswell  [30]  as  part  of  a  larger 
effort  to  parameterize  ahead  powering,  blade-rate,  fluctuating  loads,  cavitation  inception,  cavitation 
thrust  breakdown,  structural  integrity,  and  backing  performance  with  consideration  of  the  effect  of 
blade  skew.  The  cavitation  thrust  breakdown  performance  of  the  propeller  was  determined  in  a 
water  tunnel  at  a  Reynolds  number  (based  on  relative  incident  velocity  and  chord  at  70%  of  the 
span)  ranging  from  1.38  to  2.44  xlO6.  Measurements  were  taken  at  specific  advance  ratios  while 
incrementally  dropping  the  tunnel  operating  pressure.  This  process  reduced  the  test  cavitation  index 
in  a  prescribed  fashion.  Thus,  sufficient  data  were  obtained  to  complete  a  map  of  thrust  and  torque 
over  a  range  of  advance  ratios  and  cavitation  indices. 

For  the  UNCLE-M  model  cases,  the  Reynolds  number  based  on  free  stream  velocity  and 
diameter  was  held  constant  at  1.63xl06.  For  the  advance  ratios  modeled,  this  is  equivalent  to  a 
range  of  Reynolds  numbers  based  on  relative  incident  velocity  and  chord  at  70%  span  from  1 .37  to 
2.51  xlO6.  This  range  corresponds  to  advance  ratios  ranging  from  0.6  to  1.0.  The  design  advance 
ratio  was  given  by  Boswell  [30]  to  be  J=0.889.  At  advance  ratios  less  than  approximately  0.6  and 
greater  than  1.0,  the  UNCLE-M  integrations  of  the  steady  form  of  the  equations  on  the 
computational  grid  used  here  failed  to  converge.  Results  presented  here  were  integrated  based  on 
solution  of  the  steady-state  equations;  i.e.,  physical  time  derivatives  were  set  to  zero.  It  is  suspected 
that  the  lack  of  convergence  to  the  steady-state  equations  is  due  to  the  large  scale  of  unsteadiness 
inherent  in  multiphase  and  high  angle-of-attack  flows.  When  such  inherent  unsteadiness  becomes 
too  large  to  resolve  with  steady-state  physics,  it  is  typical  that  UNCLE-M  will  fail  to  converge  to  a 
steady  state  result  but  will  converge  with  application  of  the  time  dependent  formulation  [17].  In 
addition  to  requiring  unsteady  integration,  it  is  suspected  that  proper  resolution  of  the  large  scale 
structure  of  such  unsteadiness  would  necessitate  additional  and/or  better  placement  of  grid  points. 
The  grid  used  was  designed  with  a  topology  conforming  to  low  incidence  flow  without  attempting 
the  fine  resolution  of  structures  like  leading  edge  separation  eddies  and  tip  vortices.  These 
structures  were  only  captured  in  some  average  sense  and  when  the  minimum  pressure  in  such 
features  became  important,  the  overall  results  suffered.  Such  features  are  generally  not  significant 
for  breakdown  prediction  at  and  around  the  design  advance  ratio. 

The  computational  grid  used  with  UNCLE-M  to  obtain  breakdown  results  contained 
approximately  1.2  million  grid  points,  distributed  on  34  processors.  A  portion  of  the  grid  on  the 
propeller  surface  is  shown  in  Fig.  5.  Although  the  computation  was  carried  out  over  a  single  blade 
passage,  taking  advantage  of  azimuthal  periodicity,  here  the  grid  has  been  replicated  five  times  to 
give  the  appearance  of  a  complete  rotor  Results  at  design  advance  ratio  and  three  significant 
cavitation  indices  are  also  shown  in  Fig.  5.  This  picture  shows  modeled  flow  with  no  significant 
cavitation,  at  a  =1.5,  with  significant  cavitation  in  breakdown,  at  a  =1.0,  and  with  significant 
cavitation  and  in  severe  thrust  breakdown,  at  a=0.6.  Note  that  at  this  low  incidence  angle,  it 
appears  that  cavitation  is  initiated  downstream  of  the  leading  edge,  near  the  suction  peak. 

Results  are  summarized  in  Fig.  6.  As  seen  in  the  plots,  Boswell  compiled  the  torque  and  thrust 
for  the  P4381  at  a  wide  range  of  advance  ratios  and  cavitation  numbers.  It  should  be  noted  that, 
during  UNCLE-M  computations,  surface  shear  forces  were  found  to  add  approximately  5  to  1 0%  of 
the  computed  torque  and  remove  1  to  4%  from  the  computed  thrust.  These  fractions  were  higher  at 
more  lightly  loaded  conditions.  Integrated  thrust  and  torque  presented  here  include  the  shear  forces. 
In  addition  to  plotting  the  variation  in  integrated  thrust  and  torque  relative  to  significant  parameters, 
at  selected  advance  ratios  and  cavitation  numbers,  Boswell  included  pictures  and  diagrams 
depicting  the  extent  of  cavitation  over  the  propeller.  Similar  presentation  of  UNCLE-M  results  is 
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given  in  Fig.  7  (a),  for  flow  at  a  high  blade  loading  condition.  Fig.  7  (b)  includes  a  drawing  and 
photograph  from  the  original  report  by  Boswell.  Here  the  cavity  extent  predicted  by  UNCLE-M  is 
seen  to  roughly  coincide  with  the  sketch  due  to  Boswell,  only  on  the  outer  span  of  the  blade.  The 
exact  cavity  shape  is  somewhat  different.  Here,  and  at  a  higher  angles  of  attack,  leading  edge 
separation  plays  a  larger  role  in  cavity  initiation  than  at  the  design  incidence  calculations  shown  in 
Fig.  5.  Note  the  cavity  extension  into  the  tip  vortex  region.  Although  the  grid  used  here  is 
inadequate  to  resolve  tip  vortex  induced  cavitation,  it  is  clear  that  the  predicted,  attached  cavity 
persists  in  the  initial  region  of  the  (poorly  resolved)  vortex. 

In  Fig.  6  the  torque  and  thrust  breakdown  behavior  versus  advance  ratio  is  plotted  over  a  range 
of  cavitation  indices.  UNCLE-M  results  are  represented  with  open  symbols  and  data  from  the 
experiments  conducted  by  Boswell  are  represented  with  faired  curves.  Due  to  the  typical  operating 
conditions  of  a  vehicle,  of  primary  interest  is  the  behavior  near  the  design  point  and  at  reduced 
values  of  advance  ratio.  It  is  supposed  that  in  the  hypothetical  case  of  an  advance  ratio  greater  than 
design,  the  situation  would  correspond  to  windmilling  or  an  attempt  to  slow  down  or  reverse.  These 
represent  transient  events,  not  properly  addressed  here.  Therefore,  attention  is  called  to  the 
distribution  of  results  produced  by  UNCLE-M  over  the  range  of  advance  ratios  equal  to  0.9  and 
below.  For  this  range,  agreement  of  numerical  prediction  with  experimental  data  is  considered 
good.  At  every  condition  where  a  converged  solution  was  obtained,  torque  predictions  appear  to 
agree  within  5%  of  the  experimental  results.  The  agreement  with  actual  values  of  thrust  is  not  quite 
as  good  as  in  the  case  of  torque.  However,  for  more  highly  loaded  cases,  the  agreement  between 
computations  and  experiment  is  still  good. 

It  should  be  noted  that  the  actual  prediction  of  breakdown  would  be  based  on  the  prediction  of  a 
critical  advance  coefficient.  This  location  is  sometimes  based  on  some  fractional  reduction  in 
torque  or  thrust  from  the  nominal  design  value.  Typically  this  reduction  is  2-5%  of  the  nominal 
value.  It  appears  that  the  results  from  UNCLE-M,  based  on  either  integrated  thrust  or  torque,  were 
sufficiently  accurate  to  make  this  prediction.  With  proper  application  of  a  Reynolds-Averaged 
Navier-Stokes  solver,  a  high  level  of  agreement  should  be  attainable.  The  excellent  agreement 
shown  over  most  of  the  cases  for  both  thrust  and  torque  indicates  that  high  fidelity  predictions  can 
be  made,  and  with  further  refinement,  reduction  of  the  error  in  the  thrust  prediction  at  those  few 
points  where  error  was  significant  should  be  attainable. 

3.4  Multiphase,  High  Speed  Flow 

Two  applications  of  interest  are  considered  that  involve  the  modeling  of  two-phase 
compressibility  effects.  Fig.  8  shows  an  underwater  supersonic  projectile.  Both  computational 
results  and  a  corresponding  photograph  [31]  of  an  actual  test  are  included  in  the  figure.  The  three- 
dimensional,  but  axisymmetric,  computational  grid  is  comprised  of  48,049  (in  the  plane  displayed) 
grid  points,  and  second-order  spatial  accuracy  is  used.  The  flow  Mach  number  for  the  case  shown  is 
1.03  and  the  liquid  to  vapor  density  ratio  is  nominally  1000.  The  experiments  and  the  computations 
show  the  presence  of  a  bow  shock  upstream  of  the  nose.  In  addition,  because  of  the  high  velocity, 
the  cavitation  number  is  about  10  .  Consequently,  with  the  exception  of  the  nose  which  is  in 
compression,  the  flow  immediately  adjacent  to  the  body  is  completely  vaporized  as  is  the 
downstream  wake  portion. 

The  second  example  shown  in  Fig.  9,  is  the  plume  flowfield  of  an  underwater  rocket  exhaust. 
The  plume  exhaust  is  supersonic  and  is  slightly  underexpanded.  It  is  surrounded  by  a  coflowing 
secondary  subsonic  gas  stream,  which  in  turn  is  surrounded  by  a  liquid  water  free-stream  flow.  The 
nominal  liquid  to  gas  density  ratio  is  1000.  The  two-dimensional  computational  grid  is  comprised 
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of  33,153  grid  points  and  second-order  spatial  accuracy  is  used.  Fig.  9  shows  the  shock  function 
field,  which  exhibits  the  classic  expansion  pattern.  In  particular,  the  interaction  of  the  compressible 
gas  stream  with  the  incompressible  liquid  is  demonstrated  first  by  the  contraction  and  then  by  the 
expansion  of  the  gas  stream.  In  addition,  the  interface  between  the  liquid  and  gas  phases  is 
comprised  of  a  two-phase  mixture,  which  is  also  fully  supersonic  due  to  the  low  magnitude  of  the 
mixture  sound  speed. 

3.5  The  Transition  from  Reentrant  to  Twin-Vortex  Cavity  Closure 

Of  primary  concern,  is  the  nature  of  cavity  closure  behind  a  supercavitating  object.  The  mode  of 
cavity  closure  has  a  significant  impact  on  cavity  and  hence  stability.  A  frequently  encountered 
mode  of  closure  behind  large  cavities  is  the  reentrant  jet.  This  phenomenon  has  been  long  studied 
and  described  in  detail  by  Franc  [32].  Another  form  of  cavity  closure  occurs,  usually  downstream 
of  controlled  ventilated  cavities.  This  is  the  twin- vortex  type  of  cavity  closure  [2].  For  a  typical 
cavitator,  depending  on  the  Froude  Number,  the  mode  of  cavity  closure  (either  reentrant  or  twin 
vortex)  may  be  plotted  as  shown  in  Fig.  10.  Here  results  for  a  ventilated  cavity,  behind  a  conical 
cavitator,  obtained  by  Stinebring  [2],  are  shown.  In  this  figure  as  the  ventilation  rate  increases,  the 
cavity  pressure  increases  (a  decreases),  and,  to  a  certain  extent,  the  cavity  size  increases  as  well. 
Significantly,  as  the  ventilation  flow  rate  is  increased,  a  transition  point  is  found  and  the  closure 
transitions  from  a  reentrant  jet  to  a  twin-vortex  mode.  This  transition  occurs  as  a  hysteresis 
phenomenon.  Thus,  there  is  a  region  of  cavity  pressure  separating  the  two  flow  regimes  that 
appears  to  be  physically  unreachable.  In  Fig.  1 1 ,  a  photograph  of  a  cavity  in  reentrant  closure  mode 
is  shown  next  to  computational  results  at  a  dimensionless  volumetric  ventilation  flow  rate 
(Cq=q/UooDc2)  of  0.13.  The  reentrant  mode  of  closure  is  highlighted  by  both  the  overall  cavity  shape 
as  well  as  the  streamlines  shown  in  the  computational  result.  In  Fig.  12,  a  photograph  of  the  twin- 
vortex  closure  mode  behind  the  cavitator  is  shown  next  to  the  computational  result  at  a 
dimensionless  volumetric  ventilation  rate  of  0.51.  Clearly,  as  seen  in  the  photograph,  the  twin 
vortex  closure  flow  is  much  more  stable  than  the  reentrant  flow  of  Fig.  11.  The  computed  twin 
vortex  flow,  shown  in  Fig.  12,  is  also  highlighted  by  streamlines,  and  the  cavity  shape  clearly 
shows  the  twin  vorticies.  All  of  these  results  are  at  Fr=9,  the  same  as  plotted  in  Fig.  10.  Clearly,  the 
computational  results  exhibit  the  correct  mode  of  closure. 

3.6  Naturally  Cavitating  Flow  with  DES 

Fig.  13  shows  the  schematic  of  a  hydrofoil  as  offered  in  [33].  The  nominally  two-dimensional 
configuration  is  mounted  at  a  7°  angle-of-attack  and  run  at  two  different  cavitation  indices  as  well 
as  single  phase.  Though  DES  can  only  simulate  the  largest  scales  of  turbulence  when  applied  three- 
dimensionally,  in  order  to  provide  experience  and  guidance  with  the  model,  “2D  DES”  runs  were 
performed  for  the  same  two  cases.  All  of  the  simulations  (except  single  phase)  were  carried  out  in  a 
time  accurate  fashion.  Specifically,  for  all  of  the  2D  simulations,  a  timestep  of  At*  =  At/(c/U«,)  = 
0.0005  was  specified,  based  on  a  series  of  URANS  studies  performed  on  the  a=0.8  case.  A  formal 
At  study  was  not  performed,  since  as  seen  below,  even  after  60,000  timesteps,  the  simulations  are 
not  fully  stationary.  Still,  this  value  resolves  each  dominant  re-entrant  ejection  event  with  on  the 
order  of  10,000  timesteps.  For  the  2D  simulations,  10  pseudo-timesteps  are  run  for  each  physical 
timestep,  as  justified  in  [17]. 

The  inflow  velocity  was  set  to  a  constant  value  of  6m/s.  An  inlet  turbulence  intensity  of  0.01  and 
turbulence  length  scale  of  .0005*  the  channel  height  were  specified.  The  liquid  and  vapor  densities 
were  set  to  1000  and  1  kg/m3  respectively.  The  constant  liquid  and  vapor  sound  speeds  were  set  to 
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1500  and  420  m/s  respectively.  The  cavitation  number  was  defined  in  terms  of  the  average  value  of 
Pr  (see  Fig.  13),  by  accumulating  the  average  of  this  time- varying  quantity  as  the  simulations 
proceeded. 

3.6.1  2D  Results 

Five  sets  of  2D  simulations  were  performed.  A  single  phase  computation  was  carried  out  as  well 
as  two  simulations  each  at  cavitation  numbers  of  0.4  and  0.8.  The  single  phase  calculation 
converged  to  a  steady  state,  but  the  multiphase  simulations  exhibited  highly  unsteady  behavior. 
Unsteady  RANS  analysis  was  performed  for  each  of  o=0.4,  0.8.  Also,  a  “2D  DES”  simulation  was 
performed  at  each  of  these  cavitation  numbers. 

Fig.  14,  and  Table  1,  provide  details  of  the  2D  results  that  were  obtained.  The  quantities  to 
computed  are  time  averaged  lift  and  drag  coefficients,  and  various  cavity  shape  parameters.  Though 
lift  and  drag  are  unambiguously  computed,  determination  of  the  average  cavity  shape  can  be 
somewhat  discretionary.  Specifically,  in  Fig.  14,  the  time  average  cavity  shape  (defined  as  the  time 
average  contour  of  cil=0.9)  is  shown  for  each  of  the  2D  cases.  The  cavitation  detachment  location 
near  the  leading  edge  is  easily  extracted  from  this  plot.  The  maximum  cavity  length,  though 
uniquely  defined,  is  less  informative  than  when  simple  cavity  models  are  applied,  due  to  the 
complex  shape  of  the  cavity  near  its  trailing  edge.  Similar  ambiguities  apply  in  extracting  the 
maximum  cavity  thickness.  Nevertheless,  rather  than  attempt  to  smooth  or  fair  the  predicted 
average  cavity  shape  in  order  to  obtain  a  comparison  measure  more  representative  of  simpler  cavity 
models,  the  authors  have  reported  the  absolute  maxima  shown  in  Fig.  14.  The  average  total  vapor 
volume  (Vv=£avVceU  )  is  provided  in  Table  1,  for  the  four  2D  cases,  as  an  additional  indicator  of 

bubble  size.  (Reported  volumes  assume  a  spanwise  extent  of  0.1m). 

Table  1  provides  a  summary  of  the  predicted  lift  and  drag  coefficients  and  various  cavity  size 
parameters  (length/diameter,  maximum  length,  maximum  thickness,  initiation  length  to  maximum 
thickness,  all  divided  by  chord,  and  vapor  volume)  for  the  2D  runs. 


Time  average 
quantity 

Single  Phase 

ct=0.8 

a=0.4 

2D 

2D 

Cl 

0.624 

0.492 

0.171 

CD 

0.024 

0.069 

Ld/ c 

NA 

.0045 

.0077 

.0084 

Lmax/ C 

NA 

■im 

If 

1.36 

tmax/c 

NA 

0.19 

NA 

0.91 

0.73 

0.80 

BUS 

NA 

2.74 

4.38 

5.66 

7.53 

TABLE  1.  PREDICTED  CAVITATION  PARAMETERS 


3.6.2  3D  Results 

For  the  3D  DES  simulations,  the  timestep  was  doubled  to  At*=0.0010  and  the  number  of  inner 
iterates  halved  to  five.  Only  the  o=0.8  case  has  been  completed.  A  large  timestep  URANS  run  was 
first  performed  to  initialize  the  simulation.  This  was  followed  by  a  30,000  timestep  DES  run  with 
one  pseudo-timestep  per  physical  timestep.  The  “final”  DES  run  reported  here  was  run  for  30,000 
timesteps  using  5  pseudo-timesteps  per  physical  timestep.  This  required  approximately  20,000 
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processor  hours  on  a  Cray  T3E,  for  a  total  wall  clock  turnaround  time  of  about  a  week  on  100 
processors. 

The  resulting  flow  field  prediction  exhibits  a  highly  three-dimensional  unsteady  character  as 
indicted  in  Fig.  15.  There,  a  composite  of  (Xl=0.9  isosurfaces  are  presented  for  21  timesteps 
between  t*=15  and  25*.  Resolved  large  scale  features  include  spanwise,  streamwise  and  horseshoe 
vortical  structures,  a  spanwise  oscillatory  pinching  off  of  the  vapor  bubble  by  a  re-entrant  jet  in  the 
leading  edge  region,  merging  and  splitting  of  geometrically  complex  vapor  structures,  and  spanwise 
wave  structures  propagating  axially  along  the  cavity  interface  near  the  leading  edge.  Fig.  16  shows 
an  instantaneous  view  of  the  predicted  o.l=0.9  isosurface  along  with  a  number  of  particle  pathlines 
colored  by  static  pressure,  indicating  the  complex  vortical  nature  of  the  DES  flow,  including 
significant  axial  vorticity  which  vapor  structures  tend  to  track. 

3.7  Fully  Coupled  Multiphase  Flow  to  Rigid-Body-Motion 

A  three-field  (liquid,  vapor,  noncondensable  gas)  computation  of  turbulent  mixture  flow,  fully 
coupled,  via  6DOF  modeling,  to  a  supercavitating  vehicle  was  executed.  This  1.2  million  node,  48 
processor  computation,  executed  on  a  Cray  T3E  system,  was  facilitated  by  a  generous  grant  of 
Challenge  computing  resources  from  the  US  Department  of  Defense  High  Performance  Computing 
Modernization  Program  (HPCMP). 

The  modeled  vehicle  is  shown  in  Fig.  17.  Fig.  17  (a)  contains  a  schematic  drawing  and  the  basis 
for  the  notional  CFD  model.  Fig.  17  (b)  shows  a  close-up  of  the  cavitator  and  ventilation  port  area. 
This  is  from  a  steady-state  computation  of  flow  around  the  vehicle  and  illustrates  the  properly 
resolved  behavior  of  flow  around  the  cavitator  region.  The  flow  separates  and  natural  cavitation  is 
initiated  at  the  comer  of  the  flat  cavitator.  This  cavitation  is  sufficient  to  envelope  the  gas  deflectors 
which  aid  in  the  direction  of  recirculatory  ventilation  gas  exhausting  from  three  axisymmetric  ports. 
The  ventilation  gas  thus  enhances  and  stabilizes  the  cavity.  Due  to,  among  several  other  factors, 
this  recirculatory  feature,  it  is  supposed  that  a  very  low  rate  of  ventilation  flow  should  be  sufficient 
to  support  a  rather  large  cavity.  Thus  it  is  important  that  this  effect  be  captured  by  the  flow  model. 

Snapshots  of  computed  flow  field  are  shown  in  Fig.  17  (c)  Time  histories  of  vehicle  position  and 
velocity  are  shown  in  Figure  18.  For  this  computation,  the  Reynolds  Number  based  on  the  length  of 
the  body,  ReL,  is  2.3xl08,  and  the  Froude  Number  (based  on  vehicle  length)  is  13.2.  The  saturation 
vapor  pressure,  pv,  used  in  the  source  term  for  liquid  vapor  mass-transfer,  Eq  8,  is  determined  based 
on  the  cavitation  index,  a,  which  is  equal  to  0.02.  The  cavity  ventilation  rate,  Cq  is  4.4.  The  vehicle 
specific  gravity,  yg  is  1.3  with  a  uniform  density.  The  vehicle  inertial  properties  were  estimated  by 
subdividing  the  geometry  into  conical  and  cylindrical  sections  similar  to  the  description  by  Dzielski 
and  Kurdila  [34]. 

In  Fig.  17  (c),  the  vehicle  is  evident  within  the  supercavity,  and  appears  to  exhibit  a  limit  cycle 
behavior.  This  is  supported  by  the  plots  in  Fig.  18;  where  time  histories  of  vehicle  position,  velocity 
components  and  rotation  rates  are  shown.  A  limit  cycle  is  expected  for  a  vehicle  with  no  control 
surfaces  or  method  of  active  stabilization.  In  the  snapshots  showing  the  dynamic  behavior  of  the 
vehicle,  the  cavity  surface  is  illustrated  with  an  isosurface  of  liquid  volume  fraction  equal  to  0.5, 
and  the  surface  of  the  vehicle  is  colored  by  pressure.  These  snapshots  (from  dimensionless  time  1.2 
to  1 .55)  show  a  cycle  encompassing  a  complete  oscillatory  period.  This  is  clear  when  the  snapshots 
are  compared  with  the  zoomed-in  vehicle  ahead  velocity  plot  shown  in  Fig.  18  (e).  Here  a  peak  is 
shown  at  1.2  and  again  at  1.55.  It  is  at  these  peaks  that  the  high  amplitude  restoring  force  due  to 
impact  with  the  cavity  wall  occurs.  Note  the  vehicle  skips  off  of  the  lower  cavity  surface 
periodically,  and  the  impact  with  the  liquid  induces  a  high  pressure  on  the  impact  region.  This 
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impact  reduces  the  ahead  velocity  and  drives  the  vehicle  upward  slightly.  The  upward  motion, 
concurrent  with  the  reduction  in  forward  velocity,  can  be  seen  in  the  plot  of  absolute  vertical 
position.  This  sequence  explains  the  remarkable  inherent  stability  of  the  propelled,  supercavitating 
vehicle  as  configured  here  and  is  referred  to  as  planing  mode  behavior  [35]. 

In  Fig.  18  (f),  dynamic  characterization  of  the  long  term  behavior  of  the  6DOF  UNCLE-M 
model  of  the  supercavitating  vehicle  is  illustrated.  The  phase  plot  of  vehicle  relative  vertical 
velocity,  and  ahead  velocity,  is  shown  versus  time.  Although  a  simple  periodic  behavior  is  not 
expected  from  a  dynamic  model  that  has  millions  of  degrees  of  freedom,  it  is  evident  that,  based  on 
the  nearly  repeating  closed  loop  pattern,  a  quasi-periodic,  developed  limit-cycle  is  in  place.  This  is 
apparent  at  around  3.5  time  units  onward,  and  is  denoted  in  the  figure.  It  is  expected  that  more 
refined  blowing  and  propulsion  flow  rates  should  result  in  a  condition  with  much  smaller  amplitude 
limit-cycle  oscillations.  The  addition  of  control  surfaces,  such  as  cavity  piercing  fins,  should  also 
reduce  oscillations  and,  possibly,  reduce  the  vertical  drift  of  the  vehicle. 

4  SUMMARY 

A  multiphase  CFD  method  has  been  presented  and  applied  to  a  number  of  high  density  ratio 
sheet  and  supercavitating  flows.  This  presentation  begins  with  a  description  of  the  multiphase 
mixture  physical  model,  a  locally  homogeneous,  multispecies,  compressible  form,  including  mass 
transfer  modeling,  turbulence  modeling,  and  fully  coupled  flow  to  rigid  body  motion.  Several 
aspects  of  the  method  were  outlined  and  demonstrated  that  enable  convergent,  accurate  and 
efficient  simulations  of  these  flows.  These  include  a  preconditioned  differential  model  with 
favorable  eigensystem  characteristics,  regardless  of  phasic  density  ratio,  a  multiple  species 
formulation  that  separately  accounts  for  condensable  and  non-condensable  gases,  higher  order  flux 
differencing  with  limiters,  large  scale  turbulence  simulation  via  DES,  and  the  embedding  of  this 
scheme  in  a  parallel  multi-block  Navier-Stokes  platform. 

The  cases  examined,  both  validative  and  demonstrative,  show  the  capabilities  of  the 
computational  method  over  a  range  of  important  flow  conditions.  Computations  suggest  that 
UNCLE-M  has  the  ability  to  correctly  represent  the  overall  nature  of  unsteady,  complex, 
multiphase  flows.  This  in  itself  is  a  validation  of  the  approach  taken  here.  The  current  approach  has 
allowed  rendering  of  unsteady  multiphase  flows  at  Reynolds  numbers  relevant  to  engineering 
applications  in  a  modeling  method  amenable  to  complex  geometries  and  design  applications. 

The  authors  continue  to  develop  the  capabilities  of  UNCLE-M.  This  includes  the  pursuit 
relevant  validation  cases  for  complex,  unsteady,  turbulent,  high-speed,  multiphase  flows.  It  would 
seem  that  with  the  current  level  of  physical  modeling,  combined  chemistry  and  phase  change  high¬ 
speed  flows  should  be  tractable.  These  capabilities,  in  addition  to  the  already  demonstrated  abilities 
to  model  buoyancy,  ventilation,  propulsion,  and  fully  coupled  multiphase-fluid  to  rigid-body- 
motion,  are  critical  to  a  current  research  goal,  the  full  configuration,  design-quality  modeling  of  a 
controlled,  high  speed  supercavitating  vehicle  undergoing  maneuvers. 
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6  FIGURES 


Figure  1 :  Equilibrium  mixture  sound  speed  for  water  at  300K. 


Figure  2:  Computational  result.  Unsteady,  naturally  cavitating,  two-dimensional  flow.  ReL=lA*\Q5  (based  on  cavity 
length).  Modeling  of  a  two-dimensional  cavitation  tunnel.  Ref.  [28]. 

a)  Grid  with  every  4th  point  shown. 

b)  Mean  liquid  volume  fraction;  red,  a/>0.995;  blue,  a/<0.005. 

c)  RMS  fluctuating  component  of  liquid  volume  fraction;  red  indicates  a  value  of  0.5  or  greater;  blue  indicates 
negligible  fluctuating  component. 
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Figure  3:  Comparison  of  modeled,  unsteady  cavitating  flow  to  measurements  at  five  horizontal  stations  [28].  y,  vertical 
distance  from  wall,  x,  horizontal  distance  downstream  of  throat. 

a)  Mean  vapor  volume  fraction  (av). 

b)  Fluctuating  RMS  vapor  volume  fraction. 

(At  each  station,  solid  line  indicates  0  and  dashed  line  indicates  0.5.) 

c)  Mean  horizontal  velocity. 

d)  Fluctuating  RMS  horizontal  velocity. 

(Horizontal  bars  at  stations  indicate  12m/s,  the  approximate  free  stream  velocity). 
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Figure  4:  Blunt  cavitator  at  zero  angle-of-attack: 

a)  In  water  tunnel  at  0=0.35. 

b)  Model  result  from  UNCLE-M  at  a=0.4.  Isosurface  (translucent)  at  a/=0.5.  Selected  (instantaneous) 
streamlines.  Surface  of  cylinder  colored  by  a /  (red,  a/>0.995;  blue,  a,<0.005). 

c)  Experiment  to  identify  nonphysical  biasing.  Sequence  of  two  model  results  at  cr=0.35.  The  solution  in  the 
upper  left  (reflection)  of  each  snapshot  originated  from  conditions  obtained  from  the  time-dependent  solution 
in  the  lower  right.  To  obtain  the  upper  left  solution,  this  common  origin  was  transposed  to  a  mirror  image 
condition.  Isosurface  at  ar=0.5,  cylinder  colored  by  pressure  (red,  Cp>0.25,  blue  Cp<-0.3).  Strouhal  time 
indicated  \[LL*l\  from  a  reference  zero  (not  the  complete  integration  time). 
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Figure  6:  Unducted  propeller  torque  coefficient  [KQ=(Torque)/(pn2D5)]  and  thrust  coefficient  [KT=(Thrust)/(pn2D4)] 
versus  advance  ratio  over  a  range  of  cavitation  indices,  experimental  data  [30]  and  UNCLE-M  results. 


Figure  7:  Flow  over  P4381,  J=0.7,  a=3.5. 

a)  UNCLE-M  solution,  three  views  showing  surface  colored  by  pressure  and  a  gray  isosurface  of  liquid 
volume  fraction. 

b)  Diagram  and  picture  from  Boswell  [1]  indicating  extent  of  experimentally  observed  cavitation. 
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Figure  8:  Photograph  [31]  and  model  result  of  supersonic  underwater  projectile.  Model  result  contains  flow  field 
density  contours  indicating  bowshock  and  vaporous  wake. 
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Figure  9:  Cartoon  vehicle  and  3-stream,  axisymmetric  aft  flow  region.  Supersonic  gas  (air)  center  jet  (diameter=l), 
surrounded  by  subsonic  gas  (air,  at  free  stream  velocity  with  outer  diameter=2),  surrounded  by  subsonic  liquid,  water. 


Figure  10:  Experimentally  observed  modes  of  cavity  closure  for  flow  behind  a  conical  (45°)  ventilated  cavitator. 
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Figure  1 1:  a)  Photograph  of  a  ventilating  strut  mounted  45o  conical  cavitator  in  the  12”  water  tunnel 
at  the  Applied  Research  Laboratory. 

b)  Time  history  of  the  gas-on-transient  presented  as  snapshots  of  (blue)  isosurfaces  of  a /= 
0.95.  Streamlines  shown  in  red. 
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Figure  12:  Twin  vortex  cavity  closure 

a)  Photograph  of  the  aft  end  of  the  cavity  of  a  45°  conical  cavitator,  ventilating  in  the  twin  vortex  regime,  taken 
in  the  12”  water  tunnel  at  the  Applied  Research  Laboratory. 

b)  Predicted  isosurface  of  a,=  0.5  and  selected  stagnation  pressure  contours. 

c)  Farfield  pressure  contours  and  selected  cavity  streamlines  illustrating  the  twin-vortex  nature  of  the  cavity 
flow. 
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Figure  13:  Geometric  parameters  of  hydrofoil  configuration. 


v 

L. 


Figure  14:  Predicted  time  average  contours  of  aL  =  0.9  for  the  four  2D  runs  performed.  These  were  used  to  determine 
reported  cavity  geometric  parameters.  Red:  o=0.8,  URANS,  Blue:  o=0.8,  “2D  DES”,  Green:  o=0.4,  URANS,  Purple: 

a=0.4,  “2D  DES”. 
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Figure  15:  Composite  of  aL=0.9  isosurfaces  for  21  timesteps  between  t*=15  and  25*  for  DES  computation  of  o=0.8 

cavitating  hydrofoil  case. 
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Figure  16:  Isosurface  of  aL=0.9  with  particle  pathlines  colored  by  static  pressure  at  a  given  timestep  for  DES 

computation  of  0=0.8  cavitating  hydrofoil. 


dimensionless  time 
j  tUo/L=3.4 


Figure  17:  URANS-6D0F  propelled  supercavitating  vehicle: 

a)  Schematic  (basis  for  notional  model). 

b)  Close-up  of  cavitator/ventilator  illustrating  recirculatory  behavior. 

c)  Snapshots  illustrated  by  translucent  isosurface  of  constant  volume  fraction  (a^O.5)  vehicle  surface  is 
colored  by  pressure  level. 
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Figure  18:  URANS-6DOF  propelled  supercavitating  vehicle  result:  Time  history  plots,  showing 

a)  Absolute  position. 

b)  Lateral  and  vertical  absolute  positions. 

c)  Rotation  rates. 

d)  Vehicle-relative  velocities. 

e)  Ahead  velocity  spanning  a  single  period. 

f)  Phase  plot:  ahead  and  vertical  velocities  over  time. 
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